Data Visualization of WT and Nrf2 KO (-/-) BL6 PEITC or Cranberry Treated Mice 16S Microbiome Data Analysis, September 2019 Batch

Date: Fri Dec 13 12:54:45 2019
Scientist: Ran Yin
Sequencing (Waksman): Dibyendu Kumar
Statistics: Davit Sargsyan
Principal Investigator: Ah-Ng Kong

# Taxonomic Ranks:
# **K**ing **P**hillip **C**an n**O**t **F**ind **G**reen **S**ocks
# * Kingdom                
# * Phylum                    
# * Class                   
# * Order                   
# * Family     
# * Genus     
# * Species  

# options(stringsAsFactors = FALSE,
#         scipen = 999)
# # Increase mmemory size to 64 Gb----
# invisible(utils::memory.limit(65536))
# str(knitr::opts_chunk$get())
# # NOTE: the below does not work!
# knitr::opts_chunk$set(echo = FALSE, 
#                       message = FALSE,
#                       warning = FALSE,
#                       error = FALSE)
# require(knitr)
# require(kableExtra)
require(phyloseq)

Loading required package: phyloseq

# require(shiny)
require(data.table)

Loading required package: data.table
data.table 1.12.2 using 18 threads (see ?getDTthreads).  Latest news: r-datatable.com

require(ggplot2)

Loading required package: ggplot2

require(plotly)

Loading required package: plotly

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

require(DT)

Loading required package: DT

require(lmerTest)

Loading required package: lmerTest
Loading required package: lme4
Loading required package: Matrix

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

    lmer

The following object is masked from ‘package:stats’:

    step

source("source/functions_may2019.R")
# On Windows set multithread=FALSE----
mt <- TRUE

1 Introduction

C57BL/6 wild-type (WT) and Nrf-2 double-knock-out (KO -/-) mice were given 2-week microbiome stabilization process using AIN93M diet and 8 more weeks to treat with either AIN93M or AIN93M 5% PEITC diet. Fecal samples were collected weekly, immediately frozen in liquid nitrogen and stored at -80^oC. Serum, cecal, colon epithelial and whole colon tissues at week 10 were also collected for further analyses. Baseline, week 1 and 4 fecal samples were selected for 16s rRNA sequencing.

This document examines results from the WT mice samples.

We will attampt to answer the following questions:
1. Did microbiome change over time?
2. Was microbiome affected by diet?
3. Was there a difference between the KO and WT?
4. If there was a change in microbiome composition, what functional changes did it carry? What are the essential functions of the bacteria affected by the treatment and how can this be shown in vivo (metabolites, inflammation markers, etc.)?

2 Data preprocessing

2.1 Raw Data

FastQ files were downloaded from this Rutgers Box location. A total of 144 files (2 per sample, pair-ended) and a pair of undetermined reads were downloaded.

2.2 Script

This script (nrf2ubiome_dada2_sep2019_v1.Rmd) was developed using DADA2 Pipeline Tutorial (1.12) with tips and tricks from the University of Maryland Shool of Medicine Institute for Genome Sciences (IGS) Microbiome Analysis Workshop (April 8-11, 2019). The output of the DADA2 script (data_may2019/ps_sep2019.RData) is explored in this document.

3 Meta data: sample description

# Load data----
# Counts
load("data_sep2019/ps_sep2019.RData")
# Taxonomy
load("data_sep2019/taxa.RData")
taxa <- data.table(seq16s = rownames(taxa),
                   taxa)

NOTE: correction to the meta-data! (11/15/2019)

correct_samples <- fread("data_sep2019/16s metadata Sep-2019.csv")
ps_sep2019@sam_data$DSS <- correct_samples$DSS

4 Samples

ps_sep2019@sam_data$Genotype_Week <- paste(ps_sep2019@sam_data$genotype,
                                           ps_sep2019@sam_data$time,
                                           sep = "_")
ps_sep2019@sam_data$ID <- factor(paste0(ps_sep2019@sam_data$mice_num,
                                        ps_sep2019@sam_data$cage))
ps_sep2019@sam_data$TREATMENT <- paste0(ps_sep2019@sam_data$DSS,
                                        ps_sep2019@sam_data$PEITC,
                                        ps_sep2019@sam_data$cranberry)
ps_sep2019@sam_data$TREATMENT <- factor(ps_sep2019@sam_data$TREATMENT,
                                        levels = c("000",
                                                   "100",
                                                   "110",
                                                   "101"),
                                        labels = c("Naive",
                                                   "DSS",
                                                   "DSS+PEITC",
                                                   "DSS+Cranberry"))
samples <- ps_sep2019@sam_data
datatable(samples,
          options = list(pageLength = nrow(samples)))

5 Prune data

The OTUs were mapped to Bacteria (96.07%), Eukaryota (2.95%) and Archea (0.03%) kingdoms, and 75 OTUs (0.95%) undefined.

The total of 7,867 unique sequences were found. Out of those, 7,558 were mapped to bacterial genomes.

dim(ps_sep2019@otu_table@.Data)

[1]   72 7867

# Remove OTU not mapped to Bacteria
ps0 <- subset_taxa(ps_sep2019, 
                   Kingdom == "Bacteria")
dim(ps0@otu_table@.Data)

[1]   72 7558

These 7,247 belonged to 12 Phyla. 311 of the OTUs (or 4.11% of bacterial OTUs) could not be mapped to a Phyla.

t2 <- data.table(table(tax_table(ps0)[, "Phylum"],
                                  exclude = NULL))
t2$V1[is.na(t2$V1)] <- "Unknown"
setorder(t2, -N)
t2[, pct := N/sum(N)]
setorder(t2, -N)
colnames(t2) <- c("Phylum",
                  "Number of OTUs",
                  "Percent of OTUs")
datatable(t2,
          rownames = FALSE,
          caption = "Number of Bacterial OTUs by Phylum",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t2))) %>%
  formatCurrency(columns = 2,
                 currency = "",
                 mark = ",",
                 digits = 0) %>%
  formatPercentage(columns = 3,
                   digits = 2)

6 OTU table (first 10 rows)

7 Total counts per sample (i.e. sequencing depth)

t1 <- colSums(otu[, 7:ncol(otu)])
t1 <- data.table(SAMPLE_NAME = names(t1),
                 Total = t1)
t2 <- data.table(SAMPLE_NAME = rownames(samples),
                 ID = samples$ID,
                 CAGE = samples$cage,
                 TREATMENT = samples$TREATMENT,
                 Genotype = samples$genotype,
                 WEEK = samples$time)
smpl <- merge(t1,
              t2,
              by = "SAMPLE_NAME")
p1 <- ggplot(smpl,
             aes(x = SAMPLE_NAME,
                 y = Total,
                 fill = TREATMENT,
                 colour = WEEK)) +
  facet_wrap(~ Genotype,
             scale = "free_x") +
  geom_bar(stat = "identity") +
  scale_x_discrete("") +
  scale_y_continuous("Number of Reads") +
  scale_fill_discrete("Treatment") +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1)) 
ggplotly(p1)

8 Richness (Alpha diversity)

Shannon’s diversity index was calculated for each sample and ploted over time using the 7,764 from the 13 Phylum above.

p1 <- plot_richness(ps0,
                    x = "time", 
                    measures = "Shannon") +
  facet_wrap(~ genotype) +
  geom_line(aes(group = ID),
            color = "black") +
  geom_point(aes(fill = TREATMENT),
             shape = 21,
             size = 3,
             color = "black") +
  scale_x_discrete("") +
  theme(axis.text.x = element_text(angle = 30,
                                   hjust = 1,
                                   vjust = 1))
ggplotly(p = p1,
         tooltip = c("ID",
                     "value"))

p1 <- p1 + 
  scale_fill_discrete("") +
  theme(legend.position = "top")
tiff(filename = "tmp/shannon.tiff",
     height = 4,
     width = 5,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

Test if the richness changed between the baseline and the week 8.

tmp <- data.table(SAMPLE_NAME = p1$data$samples,
                  Shannon = p1$data$value)
tmp <- merge(smpl,
              tmp,
              by = "SAMPLE_NAME")
tmp$TREATMENT <- factor(tmp$TREATMENT,
                        levels = c("DSS",
                                   "Naive",
                                   "DSS+PEITC",
                                   "DSS+Cranberry"))
tmp <- droplevels(tmp[WEEK != "week1"])
m1 <- lm(Shannon  ~ WEEK*(TREATMENT + Genotype),
         # offset = Total,
         data = tmp)
summary(m1)

Call:
lm(formula = Shannon ~ WEEK * (TREATMENT + Genotype), data = tmp)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.316186 -0.091027  0.007886  0.110704  0.293230 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)                       5.94987    0.07064  84.233  < 2e-16 ***
WEEKweek8                         0.01158    0.09989   0.116   0.9084    
TREATMENTNaive                    0.14581    0.08935   1.632   0.1109    
TREATMENTDSS+PEITC               -0.03923    0.08935  -0.439   0.6631    
TREATMENTDSS+Cranberry           -0.22582    0.08935  -2.527   0.0158 *  
Genotypewidetype                 -0.54156    0.06318  -8.572 2.06e-10 ***
WEEKweek8:TREATMENTNaive          0.01181    0.12636   0.093   0.9261    
WEEKweek8:TREATMENTDSS+PEITC      0.01652    0.12636   0.131   0.8966    
WEEKweek8:TREATMENTDSS+Cranberry  0.21535    0.12636   1.704   0.0965 .  
WEEKweek8:Genotypewidetype        0.23085    0.08935   2.584   0.0137 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1548 on 38 degrees of freedom
Multiple R-squared:  0.7845,    Adjusted R-squared:  0.7335 
F-statistic: 15.37 on 9 and 38 DF,  p-value: 3.671e-10

m2 <- lmer(Shannon  ~ WEEK*(TREATMENT + Genotype) + (1 | ID),
           # offset = Total,
           data = tmp)
summary(m2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Shannon ~ WEEK * (TREATMENT + Genotype) + (1 | ID)
   Data: tmp

REML criterion at convergence: -22.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.53721 -0.47495  0.06753  0.44489  1.53874 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.01259  0.1122  
 Residual             0.01136  0.1066  
Number of obs: 48, groups:  ID, 24

Fixed effects:
                                 Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)                       5.94987    0.07064 29.76837  84.233  < 2e-16 ***
WEEKweek8                         0.01158    0.06879 19.00000   0.168  0.86814    
TREATMENTNaive                    0.14581    0.08935 29.76837   1.632  0.11322    
TREATMENTDSS+PEITC               -0.03923    0.08935 29.76837  -0.439  0.66379    
TREATMENTDSS+Cranberry           -0.22582    0.08935 29.76837  -2.527  0.01704 *  
Genotypewidetype                 -0.54156    0.06318 29.76837  -8.572 1.55e-09 ***
WEEKweek8:TREATMENTNaive          0.01181    0.08701 19.00000   0.136  0.89350    
WEEKweek8:TREATMENTDSS+PEITC      0.01652    0.08701 19.00000   0.190  0.85139    
WEEKweek8:TREATMENTDSS+Cranberry  0.21535    0.08701 19.00000   2.475  0.02291 *  
WEEKweek8:Genotypewidetype        0.23085    0.06152 19.00000   3.752  0.00135 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                     (Intr) WEEKw8 TREATMENTN TREATMENTDSS+P TREATMENTDSS+C Gntypw WEEK8:TREATMENTN
WEEKweek8            -0.487                                                                        
TREATMENTNv          -0.632  0.308                                                                 
TREATMENTDSS+P       -0.632  0.308  0.500                                                          
TREATMENTDSS+C       -0.632  0.308  0.500      0.500                                               
Gentypwdtyp          -0.447  0.218  0.000      0.000          0.000                                
WEEK8:TREATMENTN      0.308 -0.632 -0.487     -0.243         -0.243          0.000                 
WEEK8:TREATMENTDSS+P  0.308 -0.632 -0.243     -0.487         -0.243          0.000  0.500          
WEEK8:TREATMENTDSS+C  0.308 -0.632 -0.243     -0.243         -0.487          0.000  0.500          
WEEKwk8:Gnt           0.218 -0.447  0.000      0.000          0.000         -0.487  0.000          
                     WEEK8:TREATMENTDSS+P WEEK8:TREATMENTDSS+C
WEEKweek8                                                     
TREATMENTNv                                                   
TREATMENTDSS+P                                                
TREATMENTDSS+C                                                
Gentypwdtyp                                                   
WEEK8:TREATMENTN                                              
WEEK8:TREATMENTDSS+P                                          
WEEK8:TREATMENTDSS+C  0.500                                   
WEEKwk8:Gnt           0.000                0.000              

The results showed that…

9 Remove unmapped OTUs

The 311 unmapped OTUs were removed from further analysis (with 7,247 OTUs left).

ps1 <- subset_taxa(ps0, 
                   !is.na(Phylum))
dim(ps1@otu_table@.Data)

[1]   72 7247

10 Counts at Phylum level

11 Relative abundance (%) at Phylum level

Remove phyla with relative abundance of >= 1% in less than 10% of samples.

t1 <- data.table(Phylum = ra_p$Phylum,
                 `Number of Samples` = rowSums(ra_p[, 2:ncol(ra_p)] >= 0.01))
t1$`Percent Samples` <-  t1$`Number of Samples`/72
setorder(t1, -`Number of Samples`)
datatable(t1,
          rownames = FALSE,
          caption = "Taxonomic  count table",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t1))) %>%
  formatPercentage(columns = 3,
                   digits = 1)

We will remove Chlamydiae from this analysis.

[1] "Firmicutes, Bacteroidetes, Proteobacteria, Verrucomicrobia, Epsilonbacteraeota, Actinobacteria, Deferribacteres, Patescibacteria, Cyanobacteria"

XX OTUs, down from YY OTUs in the previous table.

12 Relative Abundance in Samples at Different Taxonomic Ranks

12.1 1. Class

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                          mu$Genotype,
                          sep = "_"))
p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Class,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))
tiff(filename = "tmp/wt_class_over_time.tiff",
     height = 5,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()
print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Class,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")
tiff(filename = "tmp/wt_class_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()
ggplotly(p1)

12.2 2. Order

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                                mu$Genotype,
                                sep = "_"))
p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Order,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))
tiff(filename = "tmp/wt_Order_over_time.tiff",
     height = 5,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()
print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Order,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")
tiff(filename = "tmp/wt_Order_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()
ggplotly(p1)

12.3 3. Family

NOTE: only the first 24 families had large enough counts - ploting only them.

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                                mu$Genotype,
                                sep = "_"))
mu <- droplevels(mu[levels(mu$Family)[1:24], ])

Error in `[.data.table`(mu, levels(mu$Family)[1:24], ) : 
  When i is a data.table (or character vector), the columns to join by must be specified either using 'on=' argument (see ?data.table) or by keying x (i.e. sorted, and, marked as sorted, see ?setkey). Keyed joins might have further speed benefits on very large data due to x being sorted in RAM.

13 Session Information

sessionInfo()